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Abstract. We discuss a new analytical approach to real-time evolution in 
quantum many-body systems. Our approach extends the framework of continuous 
unitary transformations such that it amounts to a novel solution method for the 
Heisenberg equations of motion for an operator. It is our purpose to illustrate the 
accuracy of this approach by studying dissipativc quantum systems on all time 
scales. In particular, we obtain results for non-equilibrium correlation functions 
for general initial conditions. We illustrate our ideas for the exactly solvable 
dissipative oscillator, and, as a non-trivial model, for the dissipative two-state 
system. 
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1. Introduction 

Strongly correlated many-body systems are challenging due to their highly non-trivial 
interplay between different energy scales. In experiment, many of their properties 
are probed by measuring their linear response to weak external driving forces. An 
additional branch of physical phenomena is currently being explored by driving 
quantum many-body systems far out of equilibrium and studying such new states of 
matter. In particular, recent progress in atomic physics has made it possible to tune 
systems of ultracold atoms at will between different interacting regimes. For example, 
a seminal experiment by Greiner et al. pQ shows collapse and revival phenomena 
in atomic Bose gases after an interaction quench, where excellent isolation from the 
environment allows the observation of such non-equilibrium behavior at remarkably 
long time scales. 

From the theoretical side, non-equilibrium quantum systems are very challenging 
since many well-established methods from equilibrium physics cannot be directly 
applied. Significant progress has been achieved in low-dimensional systems and for 
quantum impurity systems with numerical methods like NRG and DMRG: Their 
newly developed extensions called time-dependent density matrix renormalization 
group (TD-DMRG) [2] and time-dependent numerical renormalization group (TD- 
NRG)[3l EJ) were, e.g., successfully applied to interaction quenches in the Bose- 
Hubbard model [5] and to impurity spins coupled to bosonic or fermionic baths 
[1]. More recent applications of these methods also include interaction quenches in 
fermionic lattice systems in one dimension [6] and the calculation of steady state 
currents through nano-devices [7]. 

The extension of conventional analytical techniques seems to be more difficult. 
For example, the renormalization group approach, which is designed to construct 
effective low-energy Hamiltonians, is not ideally suited to understand non-equilibrium 
situations, where energy scales well above the low-energy sector can be excited. The 
situation is not much better in weakly interacting systems, where diagrammatic 
expansions are usually the first tool to analyze ground states and their elementary 
excitations, and often yield reliable results in thermal equilibrium. Nevertheless, even 
in weakly interacting systems perturbation theory is often not capable of describing 
the long-time behaviour of non-equilibrium observables, since truncation errors can 
become uncontrolled at large time scales. 

In this paper, we present a new analytical approach to the real-time evolution 
problem which merges the advantages of perturbation theory and renormalization 
group theory, and at the same time, leaves behind their shortcomings mentioned 
above. In two short previous publications, this approach has already been introduced 
and was applied to two different problems. [H [9] Here we give a somehow more 
pedagogical introduction, provide more details and discuss its general applicability to 
the field of dissipative quantum systems. In the last twenty years, quantum physics in 
a dissipative environment has played an important role in solid state physics, quantum 
optics, quantum computing, chemical and even biological systems (for a review, see 
Refs. [HI [13]). In addition, this field is especially suitable for checking the accuracy 
of our approximation scheme since many results from other approaches are known. 
However, our approach is applicable to a much wider class of problems, including also 
lattice models, see for example Ref. [9]. 
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1.1. Summary of results 

Our approach is based on an analogy to canonical perturbation theory in classical 
mechanics. We give a simple illustration of canonical perturbation theory and 
show how canonical transformations can improve perturbative expansions in real- 
time evolution problems. We show that an analogous implementation for quantum 
many-body systems is possible, based on Wegner's flow equation approach |10j . 
Independently, the same approach has been developed in the field of high energy 
physics Using this approach, we reproduce the exact solution for a quantum 

dissipative oscillator and show that efficient and precise numerical solutions of 
the analytical equations can be obtained. We also illustrate the failure of naive 
perturbation theory in this simple quantum mechanical system. 

For the spin-boson model, we obtain results that are in excellent agreement with 
known results. In the regime of weak coupling to the bath, we reproduce exact results 
for the decohcrence of a spin without tunnel splitting. For finite tunnel splitting, 
we calculate the non-equilibrium correlation function of the spin projection and also 
obtain the quantitatively correct coherent decay of the spin polarization in the weakly 
damped Ohmic regime. 

1.2. Outline 

In section[2]we give an introduction into the basic ideas of our method and its technical 
details. First we motivate our approach using an analogy from classical mechanics. 
Then we discuss the method in detail. In section G2 we apply this approach to a 
simple exactly solvable model which is useful to understand the technical details of 
our method. In section 01 we apply our method to the non-trivial spin-boson model 
to analyze the reliability of our approximation scheme. A brief summary and outlook 
concludes the paper in section [5] 

2. Real-time evolution with the flow-equation method 

2.1. Motivation: Canonical perturbation theory in classical mechanics 

There exist well-established methods to handle time-dependent perturbations in 
classical mechanics, and early attempts to handle these problems date back already 
to Newton, [M] who considered the small distortion of the moon orbit caused by the 
gravitational force of the sun. 

Much progress in this field has been motivated by the ever increasing accuracy of 
observational data for planetary motion and satellites, and the need to make accurate 
predictions based on this. The common approaches to these problems are collectively 
summarized as "canonical perturbation theory" . [14] The basic idea behind canonical 
perturbation theory can be simply illustrated by the Hamilton function of a weakly 
perturbed classical harmonic oscillator, 

tj 2 i 1 2 i 9 4 n\ 

H =2 P + 2 q A q ' ( > 

where a quartic anharmonic term with a small coupling g -C 1 perturbs the 
trajectory. We consider the initial conditions q(t = 0) = and p(t = 0) = v + §.<?^ 3 in 
the following. 
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By exploiting the smallncss of the quartic perturbation, one might be tempted to 
employ naive perturbation theory which uses a series in powers of g as an ansatz for 
the perturbed solution q(t), 

q(t)= q ^(t)+g q V(t)+0(g 2 ). (2) 
The trajectory q(°'(t) is just the solution of the unperturbed problem, 

g<°>(t) =vsm{t). (3) 
From Hamilton's equations, 



dH _ . 
dp ^ 

—q = -^ (4) 

we obtain the equation of motion 

q< 1 \t) = -qW(t)-v 3 S m 3 (t), (5) 

which has the solution 

3 v 3 

qW(t) = -v 3 sm(t) - — (sin(t)cos 2 (t) + 2sin(t) - 3icos(i)). (6) 
8 8 

This result already reveals the caveat of this approach, since the so called secular 
term 3tcos(t) yields an error growing unbounded in time. In fact, it is a well 
known general result from classical mechanics that such secular terms invalidate naive 
perturbation theory for large times. 

Naive perturbation theory can be much improved in the framework of canonical 
perturbation theory. This approach first transforms the Hamilton function to a 
suitable normal form, and after solving the equations of motion in this canonical 
basis, the normal coordinates are reexpressed through the old coordinates. In this 
manner secular terms can be avoided. To implement this idea, we first look for a 
canonical transformation of variables 



(q,p)^(Q,P) (7) 
that brings the Hamilton function to the following normal form, denoted by H , 

H = Ho + gaHl + 0(g 2 ) with tf = Ip 2 + 1q 2 . (8) 

It is easy to see that the Poisson bracket {Ho, Hq} vanishes, hence the equations of 
motion for Q and P can be solved trivially. These variables obey the initial conditions 
Q(0) = and P(0) = 0. In our example, the corresponding transformation of variables 
is 



q(t) = Q(t)-^g(3P 2 (t)Q(t) + ^Q 3 (t) 
+ 0{g% 

p{t) = P(t) + y(hP{t)Q 2 (t)+P i (t) 
+ 0(g 2 ) 



(9) 




Figure 1. We compare the different approaches to solve the equations of 
motion for the anharmonic oscillator from Eq. JTJ. The difference between the 
numerically exact solution and canonical perturbation theory according to Eq. 
mil l can hardly be noticed. Naive perturbation theory yields large errors already 
after a few oscillations, with an error that grows linear in time t. Our parameters 
are v = 4 and coupling strength g = 0.01. 



which is performed perturbatively in the parameter g. The normal form ([5]) has 
been chosen such that the equation of motion for the new variables Pit) and Q(t) can 
now be solved exactly, without producing any secular term. Using this strategy, the 
final result is 



q(t) — v sin(wt) 



3 / 5 
— gv 3 1 3 cos 2 (wi) sin(wi) + - sin 3 (wi) 

(io) 



pit) = v cos(ujt) + —gv I 5 sin (ut) cos(wt) + cos (wt) 

+ Oig 2 ) (11) 

where w = 1 + fff^o and E = p(0) 2 /2. We show a comparison of naive 
perturbation theory and canonical perturbation theory in Fig. [U which demonstrates 
the usefulness of canonical perturbation theory. 

Canonical perturbation theory yields renormalized parameters of the unperturbed 
problem, but does not lead to secular terms. This improvement becomes directly 
visible by expanding the contribution sm(cjt) in powers of g, 



sin(cjt) = sinf (1 + -gE )tj 
3 

= sin(t) + -gE tcos{t) + Oig 2 ). (12) 

This expansion generates the secular term from Eq. ([6]) occuring in naive per- 
turbation theory, making it obvious that canonical perturbation theory contains a 
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Figure 2. How to make use of the flow-equation method to implement the 
analogue of canonical perturbation theory in quantum mechanics. U denotes the 
full unitary transformation that relates the B = to the B = oo basis. 



summation over secular terms which in total yield a much improved result in compar- 
ison to naive perturbation theory. Notice that this is true in spite of the fact that both 
approaches have been expanded to the same power in the small coupling constant. 



2.2. Perturbation theory in non- equilibrium quantum mechanics 

In quantum many-body systems, the canonical way to evaluate the real time evolution 
of observables starting from some non-thermal initial state is the Keldysh technique 
|15j . which defines a contour ordered S-matrix in order to develop perturbative 
expansions for non-equilibrium Greens functions. Just as in our example of naive 
perturbation theory, secular terms can occur in any finite order of perturbation theory, 
and there is no universal solution for how to sum up these secular terms. These 
difficulties can make it very difficult or even impossible to study the transient evolution 
of observables into a steady state. 

One might wonder why an analogue of canonical perturbation theory for quantum 
many-body systems had not been developed earlier. The basic reason is the notorious 
difficulty to transform quantum-many body systems into normal form. In addition, 
the continuum of energy scales often causes often non-perturbative effects in coupling 
constants. Since the advent of renormalization group techniques, the much easier 
problem of constructing effective low-energy theories has attracted most of the 
attention. The more general problem of constructing normal forms of interacting 
Hamiltonians has been successfully treated only during the last years. 



2.3. Flow equation approach 

A general description how to transform interacting quantum many body systems into 
a non-interacting normal form has been given in 1994 by F. Wegner.[10] Let us briefly 
review the basic ideas of the flow equation approach (for more details see Ref.|17j). A 
many-body Hamiltonian H is diagonalizcd through a sequence of infinitesimal unitary 
transformations with an anti-hermitean generator r](B), 

dJ ^ 1 = [v(B),H(B)], (13) 

with H(B — 0) the initial Hamiltonian. The "canonical" generator [10] is the 
commutator of the diagonal part Hq with the interaction part Hi nt of the Hamiltonian, 
rj(B) d = [Hq(B), Hi nt (B)]. It can be formally shown that the choice of the canonical 
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generator is by construction suitable to eliminate interaction matrix elements with 
energy transfer AE — 0(1/ %/!?). Under rather general conditions, an increasingly 
energy diagonal Hamiltonian is obtained. For B = oo the Hamiltonian will be 
energy-diagonal and we denote parameters and operators in this basis by ~, e.g. 
H = H(B = oo). 

Usually, this scheme cannot be implemented exactly. The generation of higher 
and higher order interaction terms in (|13p makes it necessary to truncate the scheme in 
some order of a suitable systematic expansion parameter (usually the running coupling 
constant). Still, the infinitesimal nature of the approach makes it possible to deal with 
a continuum of energy scales and to describe non-perturbative effects. This had led 
to numerous applications of the flow equation method where one utilizes the fact that 
the Hilbert space is not truncated as opposed to conventional scaling methods. 

In Ref. [8] , these features have been exploited in order to develop an analogue of 
canonical perturbation theory in classical mechanics for quantum many-body systems. 
The general setup is described by the diagram in Fig. [2J where is some initial non- 
thermal state whose time evolution one is interested in. However, instead of following 
its full time evolution it is usually more convenient to study the real time evolution 
of a given observable A that one is interested in. This is done by transforming the 
observable into the diagonal basis in Fig. [2] (forward transformation) : 

g = [r,(B),0(B)] (14) 

with the initial condition 0(B = 0) = A, The key observation is that one can now 
solve the real time evolution with respect to the energy-diagonal H exactly, thereby 
avoiding any errors that grow proportional to time (i.e., secular terms): this yields 
A(t). Now since the initial quantum state is given in the B = basis, one undoes the 
basis change by integrating fTJJ from B = oo to B = (backward transformation) 
with the initial condition 0(B = oo) = A(t). One therefore effectively generates a 
new non-perturbative scheme for solving the Heisenberg equations of motion for an 
operator, A(t) = e tHt A(0)e~' lHt , in exact analogy to canonical perturbation theory. 



3. The dissipative harmonic oscillator 

The dissipative harmonic oscillator is a widely used toy model in the field of quantum 
optics and is also used in many other contexts. [16j It describes a quantum oscillator 
of frequency A coupled linearly to a heat bath consisting of bosonic normal modes b k . 



H = Ab% +J2"kblb k + E 

k 

+ Y,^(b + b^)(b k + b{) (15) 
k 

The operators b k fulfill canonical commutation relations 

lb k ,b{,]=S kk , (16) 

and their influence on physical properties of the quantum oscillator can be fully 
described by the spectral function 

JH=^A^( W ~ Wt ) (17) 
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In experiment, the high frequency part of J(uj) only affects the short time- 
response of the system, and thus, this part of the spectrum is usually cut off from 
J{uj). In consequence, the function J (to) is commonly approximated by a power-law 
behavior J(o>) oc lo s , u < lj c . Three different regimes of the exponent are distinguished, 
for < s < 1 the bath is called "Subohmic" , for s = 1 the bath is called "Ohmic" and 
for s > 1 "Superohmic" . 

We imagine that the system is prepared in a well-defined quantum state at some 
time to and the subsequent real time evolution of physical quantities is then defined 
by the Hamiltonian (fT5)l . The flow equation method does not restrict the class of 
possible initial states, but in this work we will consider the system-bath complex at 
time t = to be prepared in a product state 

| a) <g> |fi) (18) 

with the quantum oscillator in a coherent state 

2 * a n 

|a)=e-V y2—\ n ), <z e M (19) 

n— v 

and the heat bath in the bosonic vacuum state In such a state, the 

displacement (x) — l/y2{(b + b^)) will be finite and the effects of decoherence 
and dissipation will manifest themselves in the real time evolution of the observable 
(x(t)). The flow equation method, outlined in section [21 can solve this problem 
exactly without approximations. First, the flow equations for the Hamiltonian are 
derived. Then the real time evolution of the operators b and b' is implemented by 
their time-dependent flow equations. This will allow us to study the forward-backward 
transformation scheme in the context of an exactly solvable model. 

3.1. Diagonalization of the Hamiltonian 

As the first step of our program, the coupled form of the Hamiltonian (|15p will be 
diagonalized by infinitesimal unitary transformations. Due to the quadratic nature of 
the Hamiltonian (fT5|) . the flow equations for the Hamiltonian can be derived in closed 
form, as shown in Ref. [19] . 

Commuting the interaction part with the non-interacting boson part of (|15[) yields 
the canonical generator i](B) — [Ab^b +Yl<k UJ kb\b k ,^2 k ^k{b + b^){b k + b\)] of unitary 
transformations: 



n(B)= Y,W)A(B)(tf -b)(b k + bl) 

k 

+ J2^(B)X k (B)(b + ^)(b k -bl) (20) 
fc 

This generator leads to additional terms in the flowing Hamiltonian that violate 
the form invariance of equation (I15p . The flowing Hamiltonian preserves its initial 
form if the generator is extended by additional terms. In the following, we omit the 
explicit dependence of coefficients on the parameter B. 
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k 



J2uk\k(b + tf)(b k -bl) 



k.q 



+ £ r, hq (b k + b\){b q - b\) + Vb(b 2 - b t2 ) (21) 

k.q 

with the coefficients 

ryf 5 = -X k Af(u k ,B) 

(2) ~ 

T)l = X k uJkf{uJk,B) 

_ 2\ k \ q Au q {f ^ B) + f B)) 

U k 

1 dA 

^ =-AAlB (22) 
which are chosen such that they leave the flowing Hamiltonian form invariant. For 
this purpose the function f(ui k ,B) is still arbitrary, except for the obvious requirement 
Xk{B = oo) = 0. For our numerical evaluations of the flow equations later on, we will 
chose f(u>k,B) = — , which leads to good convergence properties in the limit 
B -> oo 

The flowing parameters of the transformed Hamiltonian are governed by the 
coupled differential equations 



dMB) -4E-fA, 



dB V 

k 



dE Q {B) 



dB ^ lk ' 

k 



dB N ' 

d\k 
~d~B 



dXk . (i) (2) 



+ 2j2vk, q X q + 2 Vb X k (23) 

The renormalization of the bath frequencies iOk will have a vanishing effect on 
time-dependent observables in the thermodynamic limit N — > oo. Therefore, the flow 
equations of these energies can be ignored for our purposes. In the limit B — > oo, the 
Hamiltonian is diagonalized and the tunneling matrix element A is renormalized to 
some value A. 

H = Ab^b + J2^kblb k (24) 

k 

In order to proceed with the real-time evolution of observables, the second step of 
the program outlined in section [5] requires the analogous transformation of the system 
operators. 
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The flow of the bosonic operator b(B) is determined by the generator (|2"Tj) and 
yields the following structure 

b(B) = /3(b + tf)+p{b-tf) 

fc 

with the flow equations 



dB 



2 Vb f3 + 2j2akVk 



Vk,qCtq 



^=-2^-2$>, fc a 9 (26) 

9 

The initial conditions are = 0) = f3(B = 0) = 1/2, a fc (_B = 0) = a k (B = 
0) = 0. During the flow towards B — > oo, the operator 6 changes its structure into a 
complicated superposition of bath operators. 

3. 2. Real-time evolution in closed form 

It is now easy to formulate an exact transformation that yields the operator b(t), t > 0, 
which is time-evolved with respect to the Hamiltonian (|15p . First, the operator b is 
trivially time evolved with respect to the Hamiltonian ([24]) as 

6(f) e^be-™ 

This operation endows the coefficients in ([25]) with trivial phase factors 



0(f) d = cos(Af) - i/3 sin(Af) 
|(f) = f 3 cos(Af) - i/3 sin(Af) 
5fe(*) = f c\k cos(cjfef) — ia/c sin(cjfef) 

Sfe(f) = a fe cos(o;fef) - ia fe sin(w fe f) (27) 



leading to the operator 



6(f) = f)(t)(b + 6 t )+/3(f)(6-6 t ) 
fc 

+ - ft £) (28) 
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Figure 3. Comparison of the exact solution Eq. I|38|l against the naive 
perturbation theory of Eq. 13411 . The secular term occuring in the second order 
perturbation expansion yields an error growing ex t, see a). In the short time limit 
t 1, naive perturbation theory becomes exact, see b). Parameters: a = -^=, 
A = 1; Ohmic bath: J(u>) = 2auiQ(ui c - ui) with a = 0.001, u) c = 10. 



The second step is to obtain the operator b(t) from the operator bit) by reverting 
the unitary transformation U used to diagonalize the Hamiltonian (I15| . formally 
represented by the relation b(t) = U^b{t)U . For this purpose, we again make an 
ansatz for the flow of the operator b(t) 

b{B, t) = 0(B, t)(b + 6 + ) + 0(B, t)(b - 6 + ) 

k 

+ a(B,t){b k -b{), (29) 

where all coefficients have both real and imaginary part, since the initial 
conditions at B = oo are given by the complex valued expressions from (|27|) . Since 
the ansatz of Eq. (f29|) is formally identical to that of (|25|). the unitary flow of b(B, t) 
can again be calculated by using the flow equations pi?)) . The operator b(t) is now 
obtained by integrating the flow equations (|26| from B = oo to B = 0, using the 
parameters from (|29|) and the initial conditions for B = oo, as given by Eq. I|27p. 
Since all transformations are unitary, the operator w(t) is readily obtained by the 
hcrmitean conjugate of Eq. (|2U1) . 

All transformations used up to now do not depend on the initial state of the 
quantum system. In calculations of time-dependent physical quantities, the operator 
b(t) can be evaluated with respect to equilibrium heat baths at finite temperature as 
well as arbitrary non-equilibrium ensembles. 

3.3. Analytical results 

1. naive perturbation theory 

X The full unitary transformation U can be expressed as a B-ordered exponential, U = 
Tb exp(J °° r](B)dB). However, this expression is only formally useful since it cannot be evaluated 
without additional approximations. 
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In analogy to classical perturbation theory as discussed in the introduction, a 
perturbative result for quantum evolution can be obtained by directly expanding the 
Heisenberg equation of motion for an operator. The exact time evolution of the 
operator x = (b + )/y2 is 

x(t) = e mt xe- iH \ (30) 

where H is given by Eq. (fl~5f . A perturbative expansion of Eq. (|30|) is typically 
performed in the interaction picture, which we define by 

xi(t) = e- lHot x(t)e lHot , (31) 

and Hq = Ab^b +^2 k u>kb\.b k is the non-interacting part of the Hamiltonian. Eq. 
(|3ip leads to the equation of motion 

^-=i[HUt),x I (t)] (32) 
which can then be expanded as 



xi(t) = x + i dn \H{ nt (ri ) , x] 
Jo 

+ i 2 [ dr x dr 2 [HUn),WUr2),xl 



+ 0(H' nt y . (33) 

Here H( nt (t) = e~ iH ° l J2 k x k{b + b r )(b k + b\)e iHot is the interaction part of the 
Hamiltonian in the interaction picture. Neglecting all contributions of 0((iJ/ Tl( ) 3 ) or 
higher in the perturbative expansion will yield a result that is correct at least to O(A^). 

Using the inital state (|18p . the first order contribution from Eq. (|3"3")l vanishes, 
since (b k ) = (bV) = 0. Likewise, all odd orders of perturbation theory vanish. 
Evaluation of the two commutators needed for the second order contribution to (x(t)} 
yields 

(x(t)) = V2a cos(At) - J d "^T^j2 

-^2 g (cOS(wt) - COS(Ai)) 



+ -sm(At)tj+0(\t). (34) 

In analogy to our example of naive perturbation theory in classical mechanics 
from section [2 a secular term occurs in the expansion from Eq. (|34[) which leads to 
large errors on long time scales (see Fig. [3]). 

2. exact solution 



We next derive a closed analytical expression for the quantity (x(t)) as a 
benchmark for a full numerical solution of the flow equations. Several other exact 
results for the dissipative quantum oscillator have been obtained by Haake and 
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Reibold.pi] According to Eq. flS5]), we have (x(t)) = y/2a/3(0,t). The coefficient 
/3(0, t) is also given by the commutator 

/J(0,t) = i[6 (t)]_, (35) 

and using the invariance of f|35[) under the unitary flow (|25[) . we obtain 

/?(0,t) = ~[6 - 6 + ,6(i)] 

= 2^(a 2 . + S*)cos(u, fe t). ( 36 ) 

k 

It is possible to evaluate this sum in closed form by defining 



= 2(^) a fe = 2(-) ~a h (37) 
what leads to the exact result for the dynamics 

POO 

(x(t)) = 2V2a / vK(w) cos(ujt)duj. (38) 
Jo 

For an Ohmic bath with J(u>) — aw6(w c -w) , the function K(u>) can be evaluated 



Sk 



as 



K(lu) = (4awA)( 167r 2 a 2 A 2 w 2 



[A 2 - c 2 + 8 a A(-^ c + ^ln(^±^))] 2 ) . (39) 
2 u>c — W 



5.^. Numerical results 



In order to illustrate the potential applications of our diagonalization scheme, we 
calculate the expectation value (x(t) by a numerical integration of the flow equations 
Eq. (|26[) . For our numerical calculation, we specify the spectral function of the bath 
as an Ohmic one with a sharp cutoff frequency u> c . 

J(lu) = au}Q(uj c -uj) (40) 

This spectral function is discretized with N = O(10 3 ) states with a constant 
energy spacing of Aui = j£. The systems of coupled differential equations (23) and 
(26) have to be solved separately for each point in time, since the initial conditions of 
Eq. (|27j) depend on time. For a finite number N of bath modes, the initial conditions 
for the bath modes have a recurrence time of T = 2irN/ui c , and it is expected that 
the error blows up rapidly at this time scale. Indeed, we could confirm this behavior 
numerically. Nevertheless, for times t < T = 0(N/uj c ), a number of about N — O(10 3 ) 
bath modes is sufficient to agree with the exact result within an error of less than 1%. 
Summing up, numerical simulations with only O(10 3 ) bath states provide excellent 
agreement with analytical solutions, with finite size effects occuring only at time 
scales of O(N). In |Appendix A[ we briefly discuss how to efficiently implement our 
transformation scheme for finite size systems. 
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Figure 4. Comparison of the time-dependent displacement (x(t) of the 
dissipative harmonic oscillator, obtained from the analytical result of Eq. J38H 
(continuous lines) and numerical integration of the flow equations (crosses). 
Different damping strengths a have been chosen, and the tunneling matrix clement 
A is almost renormalized to zero at a = 0.012, leading to a much slower oscillation 
period. The curve has been normalized to 1 by chosing a = 1/ \/2. The comparison 
demonstrates that a numerical solution of real time observables with about 1000 
bath states reproduces the analytical result already excellent with an error below 
1%. 

4. The dissipative two-level system 

The dissipative two-level system 

H = ~i ax + t 5>* + ^ + £ u " b i b k ( 4i ) 

is a fundamental model for the description of decoherence and dissipation in 
quantum systems. [13l [12] The two- state system is represented by pseudospin operators 
ai,i — x,y,z and the effect of dissipation is caused by a linear coupling to a bosonic 
bath. All bath properties can again be modeled by the spectral function J(u>) defined 
in Eq. ([T7j). 

Typically, the effect of decoherence manifests itself if the two- level system 
is prepared in an eigenstate and the coupling to the environment is switched on 
subsequently. The observable (a z (t)) describes then the tunneling dynamics of the 
initially pure quantum state. 

This problem turns out to be non-trivial and nearly all known results rely 
on approximations. The approximation scheme of the flow equation approach can 
be applied to this problem in a controlled way as shown below. This example 
demonstrates the ability of our method to treat real-time evolution in non-integrable 
models without any problems due to secular terms. 

4-1- Diagonalization of the Hamiltonian 

In order to approximately diagonalize the Hamiltonian (|4ip , we employ the following 
generator for the unitary flow: [19] 
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n(B) = ia y ]T n<?\b k + b\) +a z J2 ti ] (b k - b\) 

k k 

+ 5>M :(&* + &£)(&! - b l)-, ( 42 ) 



kl 

with B-dependent coefficients 

(v) Afc . Wfc — A ( z ) Afe w/c — A 

% = T a ^Ta' % =-T Wfe ^TA' 

^ = 2R^(^Ta + ^Ta)- ( 43 ) 

This generator has an important conceptional property that makes it different 
from the generator (f2Tj) used in section [3tl It does not leave the Hamiltonian form 
invariant during the flow, since it generates additional interactions. In this case, 
these newly generated interactions are of 0(A|), which we neglect since we consider 
the couplings X k as a small expansion parameter. For more details about the flow 
equation approach to the spin-boson model we refer the reader to Ref. |17j . Due 
to the expansion in the couplings X k , the flow equation calculation in equilibrium is 
reliable (meaning errors less than 10%) for values a < 0.2 for an Ohmic bath. 

In the above expressions, normal ordering is denoted by : ... :, which ensures 
that the truncated higher order interaction terms have vanishing expectation values 
with respect to the quantum state used for normal-ordering. In equilibrium, normal- 
ordering is therefore performed with respect to the equilibrium ground state, b k b k , =: 
b k b \' ■ +3kk'nB(k), where ns{k) is the Bose-Einstein distribution. This procedure is 
not ideal for real-time evolution of physical observables out of a non-thermal initial 
state \ipi). In order to minimize our truncation error, we define a more general normal- 
ordering procedure 



b k b k> = '■ b k b k> '■ +dkk>n B {k) + C k k' 



J k"k' ~ ■ "k u k 

Ckk' d = (il>i\b k bl,\il)i) - 5 kk >n B {k). (44) 

Here, the correct non-thermal initial state is used for normal-ordering. The flow 
equations for the Hamiltonian then read: 

^ - -2^A fc 7 7 ^ ) (l + 2n B (fc) ( 5 fefe ' +C kk .) 



= - ( w fc - A) 2 A fc + 2 XkVki 



dB 
dE 



§ The generator 1421 is again not the canonical generator in order to leave the flowing Hamiltonian 
form invariant up to higher orders in the coupling A^. It also makes use of the approximation (cr x ) = 1 
and neglects small fluctuating parts of this expectation value. 
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In the limit B — > oo, the interaction part of the Hamiltonian decays completely 
and the transformed Hamiltonian is given by 

H = -^a x + J2^blb k (46) 

k 

with a renormalized tunneling matrix element A. It has been shown in Ref. 
[T9] that A obeys the correct universal scaling behavior for Ohmic baths[12], A oc 

A(A/u c ) a ^ 1 - a \ 

4-2. Real-time evolution of operators 

The truncation scheme for the flow of the Hamiltonian (|4"Tj) can be employed in 
the same way for the transformation of the spin operators cr^, i — x,y,z. In this 
section, only the transformations of the operator a x are presented. Details of the 
transformations of the operators a y and o~ z are given in appendix |Appcndix A| 

An ansatz for the flow of the operator a x is formally given by the commutator 
[i], o~ x ]j which contains all contributions to the flow of a x which are of first order in the 
couplings Afc. For convenience, we parametrize this ansatz for the flowing operator 
o- x (B) as 

a x (B) = h(B) a x +a z J2 (xk(B) b k + Xk(B) &[) 

k 

+ a(B) + ia y J2 (MB) b k - MB) b{) (47) 

fc 

where flk and \k are related to /ifc and \k by complex conjugation. All newly 
generated terms in the differential equation d<7 ^' > = [rj(B),a x (B)] are of O(A^) in 
the coupling constants. These truncated terms are normal-ordered according to the 
convention (|44| , which finally determines the differential flow of the coupling constants 
as 

% = ~ £ fab + * fc ) + vi z) (V k + %)) 

k 

-4^4 y) C fei (x ; +X ; ) 

k,l 

d ^ = 2h^ + Y J {%Axi + Xl) + vMi - Xl)) 
I 

%=T. (^V fe + fi k ) + 4 Z \x k + X k )) (48) 

k 

with the initial conditions h(B — 0) = 1, Xk(B = 0) = fik{B = 0) = 
a(-B = 0) = 0. In the thermodynamic limit, the observable a x decays completely, 
h(B -> oo) = 0.[17l[19] 
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For the application to real-time evolution problems, it is again straightforward to 
obtain the time-evolved operator <j x (t) by evaluating <j x (t) = e lHt tj x e'~' lHt with the 
diagonal Hamiltonian H from It is easy to see that the transformed observable 

(|47[1 remains invariant under time evolution, and only its coefficients change to time- 
dependent functions 

Xk(t) = (xfe(O)cos(At) +ifi k {0)sin{At))e-^ kt 

AfcW = (Afc(0) cos(At) + i x fe (0) sin(At))e- iw **, (49) 

where a and /i remain unchanged since these contributions commute with H. 

The operator a x (t) with coefficients (|49|) can be regarded as an effective operator 
for the calculation of observables in real-time with an error of 0(Af ). Although this 
effective operator relies on a perturbative expansion of the flow equations, it is able 
to correctly describe observables on all time scales, e.g. the error remains controlled 
also for long times. This property can be understood from the analogy to canonical 
perturbation theory of classical mechanics, where the Hamiltonian function is first 
transformed to normal form. 

Using Eq. (|48[) together with the initial conditions ([4T)|). the effective operator 
a x (t) can be integrated back into the inital basis of the problem, thereby inducing a 
non-perturbative solution of the Heisenberg equation of motion for the operator a x , 
as illustrated in Fig. O Formally, this solution for the operator a x (t) is given as 

a x (t) = h(t)a x + <r z ^2 (xk(t)b k + Xkb\j 

= a(t) + io y \Pk(t)h - Pk(t)b£j (50) 

After the coefficients h(t),Xk(t),Xk(t),l^k(t),fik(t),ct(t) have been obtained, e.g. 
by a numerical solution of the flow equations, any desired correlation function of the 
operator o~ x (t) can be calculated. 

4-3. Applications 

Although it is in principle possible to obtain analytical approximations to flow equa- 
tions like (|4"5| and (|45[) . we solve the flow equations for the spin operators numerically 
in this paper. However, these numerical solutions rely on the truncation scheme we 
introduced above. In order to check the accuracy of this truncation scheme, we first 
analyze it by comparing it against an exact result. 

1. Comparison against exact results 

An exact analytical solution of the spin-boson model can be easily obtained at the 
so called pure dephasing point, where A = 0. Since er z is a conserved quantitity in this 
case, the environment can be traced out analytically. An initial state that contains 
locally a pure state \<f>) = -4=(| J.) + | f}) will become entangled with the environment 
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Figure 5. Decay of the off-diagonal matrix element Pn(t) of the reduced density 
matrix of the spin-boson model without tunnel splitting (A = 0) at T = 0. The 
dashed and dotted lines represent the exact analytical result given by Il51|l , and the 
full lines the numerical solution of the flow equations. The result demonstrates 
that the truncation scheme for the flow equations indeed is valid on all time 
scales, with an error that is systematically controlled in O(o 2 ). The bath spectral 
function has been discretized with N = 30000 states and cut off at u) c = 1. 



during the course of time. A measure for this process is given by the decay of the off- 
diagonal matrix element of the reduced density matrix, Pn(t) — (0|2V DOSOn {p(^)}|0), 
and it was shown [20] that this can be written as Pn(t) = p-n exp(— T(t)) with 

1 f°° , T / s i r u \ !- 2cos(wt) 

r(f) = -/ dwJH coth(-) (5i) 

"JO Zl LU 

Note that the quantity p-\i(t) is identical to the observable (a x (t)), which can be 
directly obtained from Eq. (|50|) . At zero temperature this observable shows a sluggish 
decay to the ground state expectation value (o~ x )gs — 0. We now compare this exact 
result with the flow equation solution at zero temperature. For the given initial state, 
the quantity p\i(t) is in the flow equation scheme given by the real function a(t) 
from Eq. ((50)) . We evaluate this quantity for Ohmic baths with a cutoff frequency 
lo c (superohmic baths are also possible, subohmic baths have not yet been treated 
successfully within our approach). The calculations show that in the regime of low 
damping strengths a < 0.1 the agreement with the exact result is excellent, see Fig. [5] 
Deviations from the exact result are on all time scales bounded by corrections of 0(a 2 ). 
The numerical results also show the small transient oscillations of the analytical result, 
which oscillate with approximately the cutoff frequency. This can be interpreted as a 
band edge effect which would vanish if a smooth cutoff function like exp(— u>/u) c ) were 
used for the bath spectrum. 

Notice that a numerical solution of the flow equations at finite temperatures (not 
shown here) is straightforward and also in good agreement with the exact result (I51|) . 

2. Decay of a polarized spin 

The formulation of the spin-boson model was originally motivated by the so-called 
quantum tunneling problem. [121 113] In this problem, the two-state system is initially 




Figure 6. Real time evolution of the spin expectation values {a x ,y,z){t) starting 
from a polarized spin in z-direction with a relaxed Ohmic bath. We used a 
damping strength a = 0.1 and the parameters N = 2000, A = 1, T = and 
u) c = 10. 



prepared in an eigenstate | |) of the pseudospin-operator, corresponding to a localized 
state in a double-well potential, which can be reduced to an effective two-state system. 
The bath, which we again take as Ohmic with a frequency cutoff u) c , is initially in 
thermal equilibrium and denoted by |bath). Thus, the initial state \ipi) of the total 
system is given by the product state 

\4>i) = | t) ® |bath). (52) 

After switching on the coupling to the dissipative environment, the time- 
dependent tunneling dynamics between the two possible states of the two-state system 
is described by the observable 

(a z (t)) = (^K(*)|^>, (53) 

Within the flow equation approach this is given by the coefficient z(t) from the 
ansatz (IA.5I) in appendix | Appendix A| We have numerically solved the flow equations 
(|A.6|) corresponding to the observable (|53p and depicted the result in Fig. [51 As 
expected, our solutions show long-time stability without secular terms analogous to 
canonical perturbation theory. For intermediate times, our curves agree well with the 
well-established NIBA approximation, [T3"l fT2j and for long-time scales t ^> A -1 the 
expectation value (a z (t)) vanishes as expected. 



3. Non-equilihrium correlation functions 



Without any additional effort we can also calculate two-time correlation functions 
based on our non-perturbative solution of the Heisenberg equations of motion for the 
spin operators. As an example we discuss the non-equilibrium correlation function 

S zz (t,t w ) = -({<T z (t + t w ),CT z (t w )} + ), (54) 

where t w is the waiting time for the first measurement after switching on the 
dynamics. Notice that in thermal equilibrium this correlation function is time 



Unitary Perturbation Theory Approach to Real-Time Evolution Problems 20 




_1 10 20 30 40 50 

At 

Figure 7. The two-time non-equilibrium correlation function S Z z(t,tw) for two 
different damping strengths a = 0.05 and a = 0.1. The full and the dashed line 
show the result for t w = 1 and the dotted lines correspond to t w = 0. Parameters: 
A = 1, ui c = 10, T = 0, N = 1000. 

translation invariant and does not depend on the waiting time t w . It is the relevant 
quantity for, e.g., neutron scattering experiments [13] . Experimental results for this 
correlator were for example obtained for hydrogen trapped by oxygen in niobium [21J, 
where protons can tunnel between two trap sites of the Nb(OH) x sample. 

For studying the non-equilibrium dynamics of the spin-spin correlation function, 
we again use the quantum state \ipi) from (f52"| with the bath at zero temperature as 
the initial state. Within the flow equation approach, S zz (t, t w ) is readily evaluated as 

S zz (t, t w ) = z(t + t w )z(t w ) + y(t + t w )y{t w ) 

+ ^2\a-k{t + t w )a k {t w ) + a k (t + t w )a k (t w )]. 

k 

(55) 

All coefficients are explained in Appendix |Appcndix A[ It turns out from 
numerical calculations (see Fig. [7] ) that the dependence on t w is only weak in the 
limit of small damping strengths and the correlations are very close to the equilibrium 
behavior. 

5. Discussion 

Our study of real-time dynamics in dissipative quantum systems was motivated by 
developing an analogous method to canonical perturbation theory as known from 
classical mechanics. 

Our results show: (i) The unitary transformation scheme of the flow-equation 
approach reproduces the well-known advantages of canonical perturbation theory from 
classical mechanics, especially the absence of secular terms in the time evolution. For 
example in the spin-boson model, physical observables can be calculated reliably on 
all time scales since perturbation theory is performed in a unitarily transformed basis, 
(ii) The identification of a suitable expansion parameter is crucial to obtain reliable 
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results. Notice that the implementation of the flow equation scheme is not restricted 
to bosonic baths or to impurity models, see e.g. Refs. [HI HI], (hi) Our method is 
particularly suitable for studying different initial states, since all transformations and 
approximations are performed on the operator level. 
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Appendix A. Additional flow equations for spin operators 

1. Derivation of Bow equations 

We briefly provide the flow equation transformations for the spin operators o~ y 
and cr 2 , that are constructed in complete analogy to the example from sectional Using 
the generator (|12"]) . the ansatze for the flowing spin components a y and o z read 

a y (B) = h^(B)a y + ia x £ x k v) (B)(b k b\) + 0(A 2 fe ) 

k 

a z (B) = h^(B)a z + io x £ x « (B)(b k + b\) + 0(X 2 k ) 

k 

(A.l) 

The flow equations for these operators are readily derived as 
£ r^x^^^cotanh^) + 8CW) 

kk' 

-2 V ^h^(B)-2j2vkiX ( i z) (A.2) 
l 



dB 

d Xl 
dB 

and 



—j^j- - "2^% Xk> (2dfe fe 'C0tanh( — )) 

kk' 

-2^h^{B)-2Y j rnkX < f ) (A.3) 
i 

Next, we solve the Heisenberg equations of motion for the operators a y and a z , 
which have the formal solution a y / z (t) = e ztH a y / z e~ ltH , with the Hamiltonian H 
given by Eq. (jlo]) . The solutions read: 
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a y (t) = U^o-y cos(-t) + h^a z sin(-i) 



a z {t) = - h^a y sin(-t) + h^o z cos(-t) 



'^ t ) + ^)a z sin(^ 

k 

|t) +^W<7,C0 B (y- 

These operators differ formally from those of Eq. (lA.ip . therefore we chose a 
different ansatz for the transformation of these operators, 

o- v /z(B, t) = o-a; ^(ia fc (£f, i)(6 fc - b\) 

k 

+ a k (B,t)(h + b l) 

+ y(B,t)a v + z(B,t)a z + 0(X 2 k ), (A.5) 

which is for both o~ v and a z identical. This ansatz yields the time-dependent flow 
equations 



^g^ = -2 Z {B,t)^ + 2Y j a l {B,t) m 

i 

dz(B,t) ^ ^ (v) W r 

d y( B ,t) -o'STn (to-sJ*) aSrJ^xWr 

— — — - 2 2_^a k (B,t)f] k -82^% a k , C kk ', 



with the initial conditions 



A' A'' 



(A.6) 



afc(t) = cos(wfct)x^ 
a fe (t) = sm(wfct)xjfe 

y(t) =ftWco8(|t) 

2(t) = U v hm(— t) (A.7) 



2 



for <T v (t) and 



&k(t) = - sin(w fc t)xfc Z) 
Sfc(t) = cos(w fc t)Xfe Z) 
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y(t) = - feW sin(-t) 

2(t) =/i (2) cos(^) (A.8) 

for <7 z (i). The effective solutions o~ y / z (t) of the respective Heisenberg equations 
of motion are now obtained by determining Eq. (|A.5|) in the case B = 0. 

2. Numerical implementation of Row equations 

In order to numerically integrate a set of O(10 3 ) x O(10 3 ) coupled differential 
equations, an adaptive step size fourth order Runge Kutta algorithm is a fast and 
accurate choice. Discretizing a bosonic or fermionic bath with equal energy spacing 
and using about N = O(10 3 ) bath states yields very accurate solutions up to time 
scales of 0(N /uo c ). In order to determine initial conditions like Eq. (|A.7[) . it is sufficient 
to integrate the flow equations up to B — 0(N 2 /u)%), where only a fraction of 0(1/N) 
of the couplings Xk has not decayed exponentially yet. In order to integrate flow 
equations with inital conditions like Eq. (|A.7|) . a standard Runge Kutta algorithm 
works not properly for large values of the parameter B, since the exponential smallncss 
of the couplings exceeds floating point precision. In our implementation, we stored 
therefore the flow of the couplings from the diagonalization of the Hamiltonian and 
supplemented it to the integration of the flow equation with time-dependent inital 
conditions. Although this procedure cannot use an adaptive stepsize in order to control 
the error during integration, it turned out to be very precise in all cases we used it. 
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